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Abstract 

applications deal with the jet breakup, leading to either of two droplet generation regimes: jet- 
ting and dripping. The cone-jet flow pattern is explored by direct axisymmetric VOF numerical 
simulation; its evolution is studied as the liquid flow-rate is increased around the jetting-dripping 
transition. As observed in other focused flows such as electrospraying cones upon steady thread 
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emission, the flow displays a strong recirculating pattern within the conical meniscus; it is shown 



to play a role on the stability of the system, being a precursor to the onset of dripping. Close 
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to the minimum liquid flow rate for steady jetting, the recirculation cell penetrates into the feed 



tube. Both the jet diameter and the size of the cell are accurately estimated by a simple theoretical 
model. In addition, the transition from jetting to dripping is numerically analyzed in detail in some 
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£> illustrative cases, and compared, to good agreement, with a set of experiments. 
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I. INTRODUCTION 

The controllable production of small flowing geometries is a crucial challenge for chemical 
engineering [Tj and bio-industry 0, [3J. Drops, bubbles, jets and recirculation cells in the 
micro-scale provide a useful platform for diverse technical applications. Here we concentrate 
on cone-jet flow patterns, looking at the streamline geometry in the jetting to dripping 
transition. Recirculating flow is shown to take place just before the transition. This may 
have been disregarded in some experimental setups, whose prime concern was the study 
of the cone and jet geometry, or the analysis of drop generation. Its study requires either 
specific flow visualization techniques or numerical simulation methods. The conditions for 
recirculation are extremely interesting, both as an indicator phenomenon associated to the 
jetting-dripping threshold, and as an attractive technological feature. 

Small droplet generation by means of co-flowing immiscible fluid streams has become 
widespread. The intrinsic smallness of the output droplets generally leads to small Reynolds 
number flows. Thus, a number of classic studies back to Taylor [4] including the recently 
blooming field of co-flowing microfluidics take the low Reynolds number assumption for 
granted. For example, a simple scheme (a straight tube surrounds a coaxial, more slen- 
der tube, two immiscible fluids being fed through each tube) has been thoroughly explored 
by Suryo and Basaran [5] using a computer simulation: a locally extensional flow sponta- 
neously develops at the tip of the elongated drop drawn by the co-flowing liquid, causing 
the seemingly continuous ejection of extremely small droplets by tip streaming under cer- 
tain parametrical combinations. This is obtained without the burden of complex geometry. 
It is hoped, therefore, that low Reynolds number co-flowing small-droplet generation may 
become a hydrodynamic standard in a near future. Different setups have been proposed 
where the flow is driven by an external straining flow. Among them, the elegant analytical 
solution by Zhang [B] points to parametrical combinations where extremely thin fluid jets, 
even down to the molecular scale, could be continuously reached. Those jets, if confirmed, 
would yield unimaginably small droplets upon breakup. 

An interesting research field is concerned with the behavior of electrified cones and drops 
of leaky dielectric fields. The problem of freely-suspended liquid droplets deforming due 
to an applied electrostatic field was examined by Haywood, Renksizbulut and Raithby [7]. 
Collins et al. [H] reported simulations and experiments supplying a comprehensive picture 



of the mechanisms of cone formation, jet emission and break-up that occur during EHD tip 
streaming from a liquid film of finite conductivity. Lac and Homsy used a boundary integral 
method to describe the axisymmetric deformation and stability of a viscous drop in a steady 
electric field [9j. In the present paper, however, only non-electrified fluids will be considered. 

In spite of their generality and tractability, low Reynolds number flows are constrained 
by the requirement that the overall flow velocity does not grow above a certain threshold to 
ensure that inertial forces remain negligible. This constraint limits the overall productivity 
of low-Re systems. Co-flowing with inertia was successfully explored, aiming at a reduction 
of the issued bubble diameter, by Oguz and ProsperettiJTU]. Subsequently, new perspec- 
tives were open by the emergence of moderate-high Reynolds number flow-focusing [H] as 
a high-productivity alternative to low Reynolds number co-flowing systems. Compared to 
other co-flowing techniques, flow focusing (FF) stands today as a mature microfluidic stan- 
dard yielding steady capillary jets or droplets whose size is well below the scale of the flow 
boundaries. As originally conceived, FF aimed at the generation [TT] of continuous steady 
micro-jets upon focusing by a co-flowing gas stream at moderate-high Reynolds number. 
Furthermore, FF was shown [12] to produce perfectly monodisperse microbubble streams 
when the co-flowing current is a liquid. A slight variation of the concept was subsequently 
introduced by Takeuchi et al. [13] to produce microbubbles. When the axisymmetric ge- 
ometry originally proposed was reduced to a planar topology [H], particularly suitable for 
microfluidics, the scientific literature production on flow focusing underwent an enormous 
boost [T51 [TH1 [T7] . In addition, axisymmetric multiple-phase FF leading to compound co-axial 
microjetsjll] has been developed by other authors to produce microcapsules in a microfluidic 
setup at relatively low [18] and moderate Reynolds number [T9]. 

The technological applications of FF were evident from its very inception. A crucial ad- 
vance resulted from the combination of massive production (high production rates of micro- 
scopic fluidic entities) and accurate tailoring. Depending on the geometry and arrangement 
of the involved fluids (a decision determining which interfaces are to be created), the nature 
of the fluids involved, gas or liquid, and the system geometry, an output including nearly 
monodisperse micro- droplets, bubbles or complex capsules can be obtained at an unprece- 
dentedly controllable rate. Surface tension becomes a paramount ally in the conformation 
of discrete (generally spherical) fluidic units. Capillary jets have long ago been observed 
to give rise to continuous drop streams at fast emission rate upon Rayleigh axisymmetric 



breakup. Here, although surface tension is negligible compared to other driving forces in the 
global scale, it becomes the main driving agent for jet instability and breakup. Obtaining 
a jet is therefore the precondition for the creation of a fluid domain with higher velocities 
and smaller dimensions at no cost in terms of control; and surface tension is free to perform 
its conformation task in this new scale. Thus, as first proposed in FF, the steady capillary 
thin jet conformed by pressure forcing by an immiscible co-flowing fluid provides favorable 
local conditions, a suitable environment for the generation of bubbles, capsules or droplets. 

A FF capillary jet is driven by three main agents: fluid inertia, viscosity and surface 
tension. Owing to the simplicity of the slender jet geometry, which asymptotically ren- 
ders all forces strictly additive in one dimension, FF can be scaled with the help of two 
dimensionless parameters: (i) the inertia to surface tension forces ratio (Weber number) 
and (ii) the viscous to surface tension ratio (Capillary number). Other classic numbers 
such as Reynolds are combinations of the former. Nevertheless, as early noticed |11] , an 
intrinsic feature of FF, namely the presence of a focusing fluid, gives rise to supplementary 
influences issuing from the correlation between the properties of the focusing and focused 
fluids. In particular, when a liquid is being focused by a gas, the gas sheath flows much 
faster than the liquid jet at the exit orifice (Fig. 1). Thus, in addition to the extensional 
viscous forces at the neck of the meniscus, transversal viscous diffusion of momentum causes 
a non-trivial axial velocity profile. Some simplifying assumptions have been adopted, yield- 
ing accurate first order solutions [TH 1201 EH]- However, they are not applicable to predict 
critical phenomena like the onset of steady jetting, or the jetting-dripping transition, as a 
function of the working parameters. These problems have been made analytically tractable 
at the expense of a drastic geometry simplification, i.e. assuming infinite jet slenderness 
[22| [231 121] • Such simplified models are predictive in a variety of situations, but FF systems 
exhibit an intrinsically three-dimensional meniscus from which the jet or the small droplets 
issue. Simultaneous modelling of the meniscus and the jet goes beyond the scope of present 
theoretical frameworks. Thus, numerical simulation or experiment are the only avenue to 
discern the physics of the fluid emission and its parametric conditions. Some further insight 
can be gained by general scaling laws. This is the approach chosen for the research presented 
here. 

Many authors have applied numerical simulation [25] to this class of problems, where 
it has supplied welcome information on the droplet dynamics of complex flows [26J. A 



significant number of studies have been proposed on microfluidic FF devices; occasional 
comparison with experimental data is provided to validate the numerical model. A liquid- 
liquid configuration for the production of microemulsions has been simulated [27] to good 
agreement with experiments |14] . Other authors have considered the microbubbling setup 
[T71 [28] . where good experimental fit is also obtained[15j. 

In this work, we make use of numerical simulation, with some experimental support, 
to study the generation of a liquid jet focused by gas in an axisymmetric FF device, at 
moderate-high Reynolds numbers. The jet diameters obtained in the simulation are in good 
agreement with our experiments and scaling laws [TT] , a fact that fully validates our hypothe- 
ses. Among other findings, we determine the flow rate at the jetting-dripping transition for 
two combinations of Reynolds and Weber numbers. We observe that either the jet or the 
cusp-like meniscus are responsible for the global instability of the system, which drives it into 
well defined dynamical cycles (global dripping). A detailed description of the flow pattern 
sheds light on the physics of the jetting-dripping transition and the peculiar appearance of 
these two regimes in co-flow problems, as opposed to faucet jetting and dripping. 

One of the key findings of the simulation is the occurrence, under favorable driving 
conditions, of a recirculation cell in the meniscus. This is in perfect analogy to recircu- 
lating meridian fluid flows observed inside Taylor cones when electrospraying liquids with 
sufficiently large values of both the viscosity and the electrical conductivity [291 EOJ- Addi- 
tionally, in experiments aiming at the production of tip streaming patterns in liquid-liquid 
two-dimensional FF (surfactant treated interface), the streamline image of fluorescent par- 
ticles seeding the flow of the internal, aqueous liquid during thread formation, was shown 
to consist of symmetric recirculation vortices [31] . In both types of motion, either driven by 
the electrical stresses acting at the cone surface or by the external focusing flow, the liquid 
flows towards the meniscus (cone) tip, along the generatrix, and away from it along the axis. 
The problem under consideration here is comparable to these other instances of recirculating 
cell, because the driving action of the gas sheath, which causes a strong tangential forcing 
at the interface, plays a similar role to either tangential electric stress[29j 30], surfactant- 
aided liquid-liquid interaction [31], or purely surfactant- driven tip-streaming. [321 [33] In this 
work, scaling arguments are developed to describe the size and occurrence of purely FF 
recirculation cells. 
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FIG. 1: Simulated boundary geometry and fluid flow domains. 



II. GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 



The axisymmetric flow-focusing device and the computational domain used are sketched 
in Fig. [I] A constant liquid flow rate Qi, flowing through a capillary tube (outer diameter 
D2 = 2i?2, inner diameter D\ = 2R\), is forced through a coaxial round orifice of diameter 
D = 2R (nozzle) located at an downstream distance H from the tube outlet. The liquid 
stream is drawn by a constant flow rate Q g of focusing gas stream discharging through 
the nozzle into a infinitely large chamber. The gas flow is assumed incompressible, in 
asymptotic consistency with the low pressure drop at the exit orifice, a condition prevailing 
when maximum droplet size monodispersity is required. Therefore, the incompressible, 
axisymmetric and unsteady Navier-Stokes equations in cylindrical (z, r, 0) coordinates are 
used to describe the time evolution of both fluids. 

Fig. 1 also shows the boundary conditions: a) at the liquid inlet, z = —zi, a Hagen- 



Poiseuille profile, Ui(r) = Vi[l— (r/Ri) 2 ] , is specified; b) at the gas inlet, z = 0, R<i < r < R 3 , 
a uniform axial flow, U g (r) = V2, is assumed. This assumption is based on the following: 
the gas Reynolds number under the conditions explored is relatively high (of the order of, 
or above 100), while most flow-focusing devices have a gas inlet length L g which is not 
much bigger than the width AR = R% — R2, so that the relevant dimensionless number for 
boundary layer development, p g U g (AR) 2 /(p g L g ) is sufficiently above unity; c) on all solid 
walls we assume no slip and no penetration u = 0; d) at the axis r = a symmetry condition 
is applied; e) the outlet discharge chamber has been modelled as a rectangular box, z = z out 
and r = r out being two open surfaces where the pressure is set to zero. This assumption is 
discussed later on. 

Note that the corresponding gas and liquid flow rates can be derived from the inlet 
velocity field: 

Q l = / 2nrU l (r)dr- Q g = / 2nrU g (r)dr. (1) 

Jo Jr 2 

Parametric studies of the dimensionless variables involved are carried out next. The 
velocity field u = (u,v) is scaled with the mean gas velocity at the nozzle V = Q g /(irR 2 ), 
while length is scaled with the nozzle radius R, time t with R/V, and pressure p with p g V 2 , p g 
being the density of the focusing gas. A single geometrical configuration is considered in this 
work, characterized by the following aspect ratios: R\/R = 0.75,i?2/-R = 1-75, R3/R = 3.5, 
H/R = 1 and L/R = 0.75. We have chosen a liquid-gas combination where 

— = 833.33, — = 55.55 (2) 

Pg Pg 

p and \x being the density and viscosity of the liquid (subindex I) and the gas (g). This 
choice is representative of the experimental jetting of air-focused water. The problem is 
governed by the following dimensionless parameters: Reynolds and Weber numbers 

*=^, (3) 

Pg 

We=e£*; (4) 

a 

a being the surface tension between the two phases. Q is defined as the flow rate quotient: 

Q = Ql/Qg- (5) 
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For a given value of Re and We we wish to analyze the formation of a steady liquid jet and 
the dependence of the flow on the quotient Q. In particular, we identify the minimum value 
of Q, Q*(Re, We), below which the liquid jet ceases to be steady and a dripping regime is 
observed in the simulation. The regime is considered to be steady (and the jet convectively 
unstable) if the liquid meniscus remains steady for a sufficiently large period of time. 

We should point out that in order to focus a jet of liquid by gas, moderate-high Reynolds 
numbers are needed. We consider in detail two different conditions for the focusing gas: 

Case 1: Re = 465.83, We = 8.137. 

Case 2: Re = 931.666, We = 32.55. 

Each case will be explored under different flow rate quotients. 

III. NUMERICAL PROCEDURE 

In order to predict the interface geometry during the time-solution, several techniques 
have been used, falling into one of three categories. These are: (i) interface tracking meth- 
ods, including a moving mesh, [M] (ii) front tracking and particle tracking schemes, [35] 
and (iii) interface capturing methods, including volume of fluid (VoF) (361 |37] and level set 
techniques. [SS] We chose a VoF method consisting of two parts: an interface reconstruction 
algorithm to approximate the interface from the set of volume fractions and a VoF transport 
algorithm to determine the volume fraction at the new time level from the velocity field and 
the reconstructed front. The basic method is robust and flexible, and is based on widely 
used VoF schemes gQl EJ 021 S3] . 

For convenience and with the aim of making our results readily reproducible for others, 
we have used the well tested commercial solver FLUENT v 6.3 (laminar unsteady) to resolve 
the discretized mass continuity, momentum conservation, and the liquid volumen fraction 
equations in the mesh depicted in Fig. 2, generated by commercial code GAMBIT in 
FLUENT v 6.2. Observe that the smallest quadrilateral elements lie between the needle 
edge and the nozzle, where the liquid meniscus is located, and in the near-axis region, where 
we expect the development of the liquid jet. The basic mesh should be sufficiently refined 
to capture, in the absence of the liquid, the strong velocity gradients experienced by the 
gas flow at the orifice region. In the grid shown in Fig. 2 the minimum cell radial and 
axial lengths are (Az) m j„ = (A r ) m j n = 0.02. Several numerical tests with smaller size mesh 
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FIG. 2: Grid of the domain under study. A denser mesh is provided in areas where the interface 
is expected to lie. To avoid numerical diffusion of the interface, the interface region is defined with 
a higher density of nodes. 



cell have shown that this level of accuracy is comfortably sufficient to describe the gas flow 
pattern for the two cases considered (Re = 465.83 and 931.666). All results presented here 
were initially computed in that mesh. In all instances where Q was very small, the results 
were recomputed in a refined mesh with quadrilateral cells in the nozzle and jet region, with 
(Az) min = (Ar) min = 0.01. Finally, only for the more difficult cases (case 2 with Q small), 
the results where recomputed in a finer grid with (Az) min = (Ar) min = 0.005. 

A factor requiring consideration is the location of the outlet boundaries. They are to 
be sufficiently remote from the nozzle to avoid the numerical reflection of pressure waves, 
since the pressure is artificially kept fixed at these boundaries during the time evolution of 
the flow. Moreover, the artificial boundary condition causes problems when the jet flows 
through the outflow boundary, and more acutely, when a string of drops flows through it. 
When a jet extends all the way to the outlet boundary of the flow domain, the pressure 
within the jet and the surrounding fluid cannot be equal owing to surface tension. However, 



for our particular problem, as will be shown later, these undesired effects are confined 
to a length below two diameters upstream of the boundary. Therefore, we have chosen 
simulations with a sufficiently large external chamber, z out = 10 and r out = 3.5, to minimize 
artificial boundary effects in the results obtained. On the other hand, the z position where 
the inlet boundary conditions for the liquid are imposed, has been located sufficiently far 
away from the needle edge, at z = Zi — —3. This choice has been made, as will be shown 
later, because a liquid recirculation cell intrudes upon the capillary tube when Q decreases. 
Therefore, in order to impose a Hagen-Poiseuille profile for the liquid velocity as a well 
posed inlet boundary condition, this boundary should be set sufficiently far upstream from 
the recirculating region. 

Tracking the interface between the phases is accomplished by solving a continuity equation 
for the volume fraction of one of the phases using an explicit time-marching scheme. The rest 
of the equations are solved implicitly. The time steps selected were fixed and sufficiently 
small to ensure that the global Courant number based on the mesh cell size, the mean 
velocity in the cell and the time step was always smaller than one. Regarding the spatial 
discretization of the equations, the third-order modified MUSCL scheme |B] is used to obtain 
the face fluxes whenever a cell is completely immersed in a single phase. When the cell is near 
the interface, the CICSAM algorithm is used[45j. The pressure corrections are computed 
with the body forces weighted scheme and the pressure-velocity coupling in segregated solver 
is treated with the PISO method |46j. All under-relaxation factors are set to one to avoid 
any numerical masking of fade-out effects in our physical problem. 

IV. NUMERICAL RESULTS. DISCUSSION 

A fruitful interpretation of the results obtained needs to be situated in the frame of the 
literature on the dripping faucet. The book by Shaw |Hj gave rise to a rich and insightful se- 
ries of studies, among them major contributions by Fuchikami et al. [48J, Ambravaneswaran 
and co-workers [H2 [50], [51] , and Coullet, Mahadevan and Riera. [52] To discriminate between 
the jetting and dripping modes, it is helpful to make use of the categories introduced by 
Ambravaneswaran and co-workers: 



1. the dimensionless limiting length L^/Ri from the capillary edge to the extremity of 
the first drop at detachment. 
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2. the ratio of the distance L s between the centers of mass of the drop that is about to 
form and the previously formed drop, and L^ 

3. the ratio of the volume of the drop that is about to form, Vd, to that of the drop that 
is attached from the capillary, V p 

When undergoing the transition from dripping to jetting, the first parameter undergoes 
a sudden increase, while the two other ones experience an opposite trend. In general, the 
dripping mode is characterized by bulky drops, relatively distant from each other, whose 
diameter is considerably larger than the jet diameter from which they detach. 

The usual categories applicable to faucet dripping need some adaptation before being 
used in a co-flow problem, where there is considerable stretching of the cone-jet and droplet 
train by the coaxial current: the drops are deformed, their radial extent is limited, and 
they may undergo secondary breakup (particularly so in the dripping regime, whose bulky 
drops are more vulnerable to shear) after detachment from the filament. Therefore, under co- 
flowing, the classical aspect of the jetting and dripping regimes is modified, and the transition 
between them is not sharp: such features are confirmed by experiment, as explained below. 
This is the reason why the behavior of the meniscus can be used as a further indicator of 
the jetting mode. Dripping leads to a pulsating meniscus, each detached drop giving rise 
to recoil and oscillating; while, in jetting, the detachment of the drops does not cause any 



fluctuation of the meniscus and jet (see Fig. 19). In a full dripping regime, these pulses are 



perfectly regular (see Fig. 14); in most cases a slender unsteady liquid ligament detaches 
from the meniscus and breaks up into droplets of heterogeneous size [53J. However, at 
the onset of dripping (a situation which will be labelled "incipient dripping"), completely 
irregular fluctuations of the meniscus are observed |5Uj . 

Accordingly we begin by studying the formation of a steady (convectively unstable) liquid 
jet in the FF device. Initially, the capillary needle is filled with liquid up to z = while the 
rest of the domain is filled with gas. We start the simulation from rest (u — v — 0) in the 
whole domain except at the inlet sections, where velocity profiles are prescribed. Fig. [31(a)- 
[3]^h) shows the formation of a steady liquid jet for case 1 and Q = 0.004, going through the 
stages of interface entry, meniscus growth and jet consolidation. The shape of the liquid-gas 
interface is computed in the figure as the iso- level of the liquid volumen fraction a = 0.999, 
obtained with the VOF method. Given that the flow might be unsteady, we consider that 
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the meniscus-jet has reached a steady condition whenever two conditions hold: 

1. The angle between the liquid meniscus and the radial coordinate at the capillary 
needle, 9(t), has reached a constant value in time, 9{t) = 9 . 

2. Both the jet diameter at the nozzle inlet, di n (t), and at the nozzle exit, d ou t(t), should 
reach a steady regime or a stable oscillating regime around a mean value; these quan- 
tities are of course to stay above zero. This amounts to excluding jet breakup in 
the nozzle region, a feature associated with a non-slender jet and possible dripping 
behavior (unsteady meniscus-jet). 

Here, d in (t) and d out (t) are computed at each time step, by integrating radially the liquid 
volume fraction, a, at the nozzle inlet, z — 2, and at the exit, z = 2.75: 



(/,„(/-) = 2\ / 1 I a(t,z = 2,r)rdr; d out (t) = 2W2 f a(t, z = 2.75, r)rdr. (6) 

For sufficiently large Q, as illustrated in Fig. |3j both d in and d out evolve towards a steady 
value. 

However, oscillations of these two quantities are observed when Q is reduced. For exam- 
ple, Fig. [4] shows the time oscillation of d in and d out for case 1 and Q = 0.0006 after allowing 
a steady jet to develop. It can be observed that the jet diameter at the nozzle exit is smaller 
than at the inlet; mass conservation arguments imply the inlet velocity to be smaller than 
the outlet velocity (in inverse proportion to the diameter squared). This explains why the 
oscillation frequencies of the jet diameter are shorter at the outlet. Although the oscillation 
of the jet in the nozzle region may play an important role in the dynamics of the droplets 
generated upon jet breakup, our main concern here is to characterize the jet diameter, the 
angle of the meniscus at the attachment, 9, and the flow structure inside the meniscus as a 
function of Q. Since the flow is unsteady, we will use a mean value of the jet diameter at 
the nozzle inlet and outlet as defined by: 

i pT+U _ i pT+ti 

d in = — d in (t)dt, d^ = 7p I dout(t),dt (7) 

where U is a time position once a steady jet has developed and T is a time period long 
enough to ensure a significative mean value. For example, selecting t, = and T = 500 
leads to d in = 0.2456 and dout = 0.1527 (conditions as in Fig. |4j). 
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FIG. 3: A sequence of snapshots from the simulated growth of an eventually steady jet (case 1, 
Q = 0.004): (a) the interface arrives to the needle edge; (b)-(c) the meniscus grows in the nozzle 
region; (d)-(f) a jet begins to issue from the nozzle; (g)-(h) the meniscus-jet system is steady. 

The procedure is the same for the two cases under consideration. The simulation is 
started from rest with a value of Q sufficiently high to obtain a steady jet. Then, Q is 
reduced and the solution is monitored in time until a new steady jet is obtained. Fig. [5] 
shows the stabilized liquid-gas interface for case 1 and different Q. It should be pointed out 
that Q = 0.0004 is the smallest flow rate compatible with a steady jet for case 1. Therefore, 
it can be identified as the minimum flow rate Q* for steady jetting: Q* = 0.0004 for case 1. 

Fig. [6] shows the interface for case 2 and different Q values once a steady regime is 
reached. The smallest jetting flow rate here is Q* = 0.0001, four times smaller than in case 
1. Accordingly, the smallest jet diameters are obtained for case 2. The jet diameter evolution 
is shown in Fig. [7], where the mean steady values di n and di n are plotted as a function of 
Q for (a) case 1 and (b) case 2. To complete the picture, Fig. M shows the dependence of 
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FIG. 4: Time dependence of the jet diameter, di n (dashed line, left ordinates) and d ou t (continuous 
line, right ordinates) for case 1 and Q = 0.0006. Note (i) the smallness of the oscillation amplitudes 
and (ii) the higher frequency in the variation of d ou t as compared to di n . 

the meniscus angle 9 with Q for the two cases. In both examples, 9 becomes smaller as Q 
decreases (smaller flow rate quotient implies stronger focusing action). Just before dripping, 
as the liquid flow rate is reduced, the angle appears to become independent from Q: the 
interface geometry becomes invariant (local hydrostatic balance). The smallest angles are 
obtained in case 2. This is to be expected since the normal pressure forces produced by the 
gas stream, which cause the focusing flow, are larger for that case. 

Analyzing in more detail the structure of the flow inside the liquid meniscus in jetting 
mode, in the lower-Q range, a meniscus recirculation cell is observed, in analogy with other 
co-flowing systems [3 [31] and Taylor cones [29] . Fig. M shows instantaneous streamlines 
for case 1 and different values of Q. The recirculation increases when Q decreases, the cell 



penetrating into the capillary needle. Fig. 10 depicts instantaneous streamlines for case 2 



and four different values of Q. Again, a recirculation region appears before the meniscus jet 
system ceases to be steady. The size of the recirculation can be calculated by finding the 
two ^-positions where the velocity at the axis becomes zero. 
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FIG. 5: Shape of the liquid-gas interface as a function of Q (case 1) 



Fig. [TT] shows the velocity at the axis, v aX { s , as a function of z for five different values 



of Q: (a) case 1 and (b) case 2. It is worth observing that v ax i S is roughly uniform inside 
the capillary needle, its value being given by the Hagen-Poiseuille expression; in the nozzle 
region it increases owing to the focusing effect of the gas stream, which creates the issuing 
jet. A region where u ax i s decreases is located in the meniscus region, between the capillary 
and the nozzle; note that when Q decreases, a local minimum of u ax i S is observed at a given 
position z = z m i n in the meniscus region. If Q is sufficiently small, v ax i S becomes negative 
near the local minimum in a region delimited by the two z positions, z\ and z 2 , where 
Vaxis — 0. Therefore, the size of the recirculation region, S r , observed in Figs. [TO] and [TT 



can be computed as S r /R — s r — z-i — Z\. There is a threshold value of Q, Q r , below which 
a recirculation pattern is observed. At the threshold flow rate v ax is — at z\ = z m i n = z 2 
and v ax i S > elsewhere. 



Fig. 12 shows s r as function of Q for (a) case 1 and (b) case 2. Looking back at Fig. pi 
note that the size of the recirculation cell increases as Q decreases. In situations of incipient 
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FIG. 6: Shape of the liquid-gas interface as a function of Q (case 2) 





FIG. 7: Mean jet diameter at the inlet (dashed line) and at the exit (solid line) of the nozzle orifice 
versus the flow rate quotient Q. (a) Case 1, (b) Case 2. 
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FIG. 8: Angle 6 (degrees) versus flow rate quotient Q. (a) Case 1, (b) Case 2. 
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FIG. 9: Instantaneous streamlines at four different flow rates quotients Q for case 1: (a) Q = 0.004, 
(b) Q = 0.0008, (c) Q = 0.0006 and (d) Q = 0.0004. 



recirculation (Q smaller than but similar to Q r ) this growth appears to be linear, as derived 
later from dimensional arguments. In Fig. [12J the discrete points V have been obtained 
directly from the simulations. The dashed lines are linear interpolations computed in the 
recirculating regime, s r > 0. The linear interpolation is not only in good agreement with 
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FIG. 10: Instantaneous streamlines with four different flow rate quotients Q for case 2: (a) Q 
0.0012, (b) Q = 0.0004, (c) Q = 0.0002 and (d) Q = 0.0001. 
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FIG. 11: Velocity at the axis versus z for several flow rate quotients Q: (a) Case 1 and (b) Case 2. 
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FIG. 12: Size of the recirculation s r as a function of Q: (a) Case 1 and (b) Case 2. 

the data but also provides a reliable approximation to compute Q r . The estimations are 
Q T = 0.001453 for case 1 and Q r = 0.000708 for case 2. According to the above, the linear 
expression relating the size of the recirculation region and the flow rate Q is: 



s r ~A{Q r -Q). 



(8) 



This means that s r is proportional to a back flow rate QbIQq = Qb = {Qr — Q) for a 
given set of fluid properties, geometry, and gas flow Reynolds number. The relative location 
of the jetting threshold Q* and the recirculation threshold Q r is, in both cases, Q* < Q r . 
Therefore, recirculation can be taken as a dripping-precursor: as far as can be gathered from 
our simulation, it always precedes global instability of the meniscus-jet system. 

Finally, some results are presented with flow rate values below the jetting threshold. For 
Q < Q*, with Q close to Q*, our simulations show the flow to exhibit different behaviors in 
a sequence: a period where a thin jet breaks up in the nozzle region alternates with other 
periods where a thin jet breaks up downstream of the nozzle. The irregular time behavior 



of the flow for Q < Q*, but Q close to Q* (incipient dripping) can be observed in Fig. 13 



where di (a) and d out (b) are shown as a function of time for case 1 and Q = 0.000322. For 
Q < Q* but Q sufficiently different from Q*, the flow behavior becomes more regular and 
periodic with a unique dripping frequency. 

As anticipated at the first part of this section, our results, while belonging to the dripping- 
regime, are untypical in that the co-flowing current gives rise to axial stretching of the 
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FIG. 13: Time evolution of (a) di n and (b) d ou t in a irregular dripping regime for case 1 and 
Q = 0.000322. 





FIG. 14: A sequence of the dripping regime for case 1 and Q = 0.000241 involving a complete 
cycle. 
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FIG. 15: Time evolution of (a) di n and (b) di n in a dripping regime for case 1 and Q = 0.000241. 

jet and drops, so that unusual breakup geometries result. The radial extension of the 
drops is limited, and the axial stretching gives rise to secondary breakup, immediately after 
detachment. The pattern observed in Fig. [14] points to a dripping regime: it is perfectly 
periodic (transient jetting can therefore be excluded) and each period is associated with the 
filling up of a drop, its breakup from a thinning filament, and the recoil of this filament. Fig. 
14 shows a complete time sequence of a dripping process (case 1, Q — 0.000241). Subfigures 
(a)-(e) show the growth of the meniscus and the formation of a jet issuing from the nozzle; 
subfigures (g)-(h) show the jet breakup into droplets of different sizes and the meniscus 



recoil. This is a periodic sequence, the period being T ~ 210 for each cycle. Fig. 15 plots 
di (a) and d out (b) as a function of time for this case. Initially, a liquid meniscus is growing 
with no jet production, and d{ n = d out = 0. Subsequently, a liquid jet issues and di n and 
dout become positive. Both diameters reach a maximum at a certain time and then start to 
decrease. Finally, the jet breaks into droplets, di n and d out are set to zero and the process 
begins anew. In spite of the observed differences, the dripping process in this case is quite 
similar to regular faucet dripping, the time period being mainly imposed by the filling of 
the meniscus until reaching a critical volume. 
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FIG. 16: A sequence of the dripping regime for case 2 and Q = 0.00004 involving a complete cycle. 



A similar situation is observed in case 2. Fig. 16 shows a complete time sequence of a 
regular dripping process with Q = 0.00004. Subfigures (a)-(d) show the meniscus growth 



and the emission of a jet, much thicker and longer than observed in Fig. 14 Again, this 
sequence is time periodic with a period T ~ 500 for each cycle (see di (a) and d out (b) in 



Fig. 17 as a function of time). 



A. Influence of the BCs and the spatial and temporal resolution on the numerical 
results 

The numerical problem addressed is quite complex: it involves a high speed stream of 
gas discharging through a nozzle into a infinity large chamber plus a meniscus-liquid jet 
which may break into droplets within the finite numerical domain. This complexity leads to 
different time and spacial scales associated to a plurality of interacting physical phenomena 
(jet breakup due to capillary and Kelvin-Helmholtz instabilities, mixing layer instabilities 
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FIG. 17: Time evolution of (a) di n and (b) di n in a dripping regime for case 2 and Q = 0.00004 

in the main gas stream). Therefore, an accurate analysis of the jet breakup is difficult to 
achieve. On account of it, though our VOF method is fully reliable in qualitative terms as 
a predictor of jet breakup and drop formation, we have devoted this subsection to check 
that our numerical results (meniscus-jet shape as a function of Q for different setups) are 
independent of the selected BCs and numerical meshes. 

As indicated above, the most problematic simulation choice is setting p = at the 
outlet boundaries, since any jet or a drop crossing the boundary is influenced by the strong 
and artificial restriction that the pressure remains fixed. Our choice is a simplification 
(p = 0) which takes advantage of the essentially parabolic character of the equations. A 
second option has been explored, the so-called outflow conditions (assuming uniformity, i.e. 
Neumann type), but they give rise to a false constraint on the flow pattern, because they 
imply that the gas flow is coaxial. A minisymposium held in 1994 on the open boundary 
condition problem in incompressible flow, by Sani and Gresho [54J led to the concluding 
remark: "We have made some attempts at shedding more light on the difficult and unresolved 
area of seeking good OBCs for incompressible flow simulations. It has been an exercise in 
frustration and we are not thrilled with the results obtained" . 
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FIG. 18: The effect of the p = boundary condition on the liquid-gas interface comparing to 
downstream boundary locations: z = 6 (thicker line) and z = 10 (thinner line): (a) data as in case 
1 and Q = 0.004; (b) data as in case 1 and Q = 0.0008. 

There is an evident inaccuracy involved in our p = choice: the pressure jump associated 
to an interface will lead to high local pressure inside the liquid jet or droplets. However, 
this assumption can be reconciled with our aim, which is not a study of the breakup process 
and its transient geometries. We are addressing a wider scale: the cone-jet flow pattern, and 
the general drop generation regime. To show that the distortion caused by this artificial BC 
is local and does not modify the global behavior at the cone-jet region, some exploration 
as been carried out. It can be shown that setting the external boundary sufficiently far 
downstream from the nozzle region, at z out = 10, the meniscus-jet is not affected by the 
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FIG. 19: The effect on the liquid-gas interface of the spatial resolution: (a) a general view of case 2, 
Q = 0.0001 where the grey line is the solution computed in a mesh where (A 2 ) m j n = (A r ) m j n = 0.01 
and the black line is the solution obtained in a mesh with (A z )mi n = (A r ) m j ra = 0.005; (b) contours 
of the vorticity field for the same case with the best spatial resolution. 



boundary condition. To show this, we have considered the worst scenario: we choose large 



liquid flow rates Q and weak gas flow (case 1). Fig. 18 shows the stabilized liquid-gas 
interface for case 1 and two different values of Q, computed in the original domain and in 
a shorter one. In Fig. fl8ta), the jet does not break up within any of the two numerical 
domains and the jet and meniscus interface in the nozzle region is evidently not affected by 
the artificial p = boundary condition. The influence of the artificial BC is confined to a 
few diameters upstream of the downstream boundary. In the case considered in Fig. [i~8J(b), 
the jet is breaking up into drops within the large domain. Even in this case, the meniscus 
and jet interface in the nozzle region is not affected by the artificial boundary condition. 

Let us now show the consistency of the model by comparing the results in two different 
meshes. In this case again, we have considered the worst scenario, by selecting smaller values 
of the liquid flow rate (Q small) and a large gas flow (case 2), since thinner jets are obtained 



25 




120 




FIG. 20: a) Comparison of the time evolution of the vorticity magnitude at the point {z = 9.5, 
r = 0.995) for case 2 and Q = 0.0001 computed with At = 0.014 (solid line) and with At = 0.028 
(dashed line); b) shows the frequency content of the two signals. 



in these cases. Fig. 19 (a) shows a instantaneous picture of a steady (convectively unstable) 



liquid jet breaking up into drops (case 1, Q — 0.0001) computed in two different meshes. 
Observe that although we are comparing the liquid interface at two different times and 
with different spatial resolution, the shape of the meniscus of the liquid in the nozzle zone 
coincides. The main difference is that the liquid jet is slightly longer in the finer mesh. As 
mentioned, an accurate description of the jet breakup is not the objective of this paper. It 
requires specific analytic tools in order to capture the diverse physical phenomena involved. 



This can be illustrated by Fig. 19 (a), showing instantaneous contours of the vorticity 



field for this case, computed with the best spatial resolution. It can be seen that, owing to the 
large density and velocity difference between the liquid and the gas, the vorticity is large in 
the liquid-gas boundary layer which develops at the liquid-jet meniscus. When the jet breaks 
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up, the behavior of the flow around the drops is similar to high Reynolds number flow around 
a round object (since the gas is travelling faster that the liquid drops). Therefore, boundary 
layer separation is to be expected at the surface of the drops, and a wake will originate 
downstream of the drop. In addition to this, a mixing layer is developing downstream of 
the nozzle edge owing to the difference between the velocity of the gas flowing through the 
nozzle and the stagnant atmosphere surrounding the outlet. This mixing layer yields to 



the development of vortices as shown in the figure. Fig. 20 a) shows the time evolution of 
the vorticity magnitude at a point of this mixing layer computed with two different time 
steps. Strong fluctuations of the vorticity magnitude are apparent. The flow evolution shows 
some sensitivity to the resolution, particularly as the simulation time increases. However, 
both time resolutions are reliable to predict the characteristic frequencies of the problem at 
the observation point. Fig. |20|(b) show the frequency content of the vorticity as obtained 
by applying the Fast Fourier transform (FFT) to the two signals. The main frequency 
u>i ~ 0.46 is related to the passage of the vortices generated in the gas mixing layer. There 
is a secondary characteristic frequency, uo 2 ~ 0.56, associated to the interaction between the 
mixing layer and the vorticity wake of the drops. This interpretation was strengthened by 
recomputing the flow without the liquid jet: the frequency content of the vorticity signal at 
the same observation point only showed the main frequency peak uj\ ~ 0.46 . 

V. COMPARISON WITH ANALYTICAL MODELS AND SCALING LAWS 



The first predictive model for the jet diameter dj at the orifice exit [TT] assumes that 
viscous and capillary effects are small enough compared to liquid inertia. This demands 
large enough Reynolds and Weber numbers of the liquid jet, in reasonable agreement with 
most experimental conditions (common solvents including water, down to the micron scale). 
In this limit, the overall pressure difference AP = AP(Q g ) (pressure difference between the 
gas inlet and the gas outlet) imposed in the downstream direction (i.e. through the orifice), 
transmitted to the liquid stream by normal surface stresses, is converted into kinetic energy, 
so that 

A H"^S (9) 
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FIG. 21: Jet size measured at the entrance of the nozzle using d{ n (□ correspond to case 1 and 
* to case 2) and at the exit using d ou t (o correspond to case 1 and x to case 2) compared to the 
theoretical prediction (continuous line). 



which readily gives 



dj 



i/i 



Q 



1/2 



(10) 



Furthermore, the jet is assumed sufficiently small compared to the orifice diameter D such 
that no only it does not touch the orifice borders, but also the boundary layers of the 
focusing fluid (gas) at the orifice and at the jet's surface are sufficiently small compared to 
the corona defined between the jet and the orifice. This is why D does not enter expression 



(10). Neither does D\ have any direct influence on the jet diameter; only as a parameter 



determining the liquid flow rate Qi. 

Interestingly, if viscous effects and surface tension are neglected, and we assume dj « 
D, the only operating parameters left in the analysis are {pi,AP, Qi}\ using these three 



parameters, a scaling law identical to 10 -regardless the constants- follows from dimensional 



reasoning. Figure 21 illustrates the accuracy of this first simple prediction. However, that 
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expression does not provide information on how small the neglected effects are. Would it be 
possible to quantify both the dependence of the jet diameter on the three main parameters 
and the relative magnitude of each one of the neglected effects? The answer is yes, by 
retaining either a, //j, or D in the dimensional analysis, and using {AP, a, pi}, {AP, Hi, pi} 
or {AP,D,pi} as independent parameters, respectively. 

In physical terms, the relative effect of surface tension may be determined by observing 
that the liquid Weber number must be of order unity (Wei = Pi&Qf/('X 2 0'd%)) for a given 
pressure AP, from which one obtains the limiting diameter d and flow rate Q . The result is 

1 /9 

d = o J AP and Q = (cr 4 p z -1 AP~ 3 ) , so that dimensional analysis together with equation 



(10) readily yields 

d 3 /d = k d {Q l /Q ) 1 ' 2 (11) 

where kd = (8/7r 2 ) 1//4 . This expression provides a first-order approximation to the jet diam- 
eter (asymptotically true for Re —* oo) as far as dj » d (negligible surface tension). The 
ratio Qi/Qo spans the whole domain from jetting to dripping -where dj becomes comparable 
to d -. Jet diameters and flow rates comparable to d and Q , respectively, lead to abso- 
lute instability, where the characteristic velocity of upstream capillary waves 0(apY 1 dJ 1 ) 1 ^ 2 
-a product of surface tension- becomes of the order of the downstream convective velocity, 



Qi/d 2 . Besides, Eq. (11) is explicitly independent of the orifice diameter P, an illustration 
of the jetting regime independently of its forcing geometry. 

Similarly, viscous effects can be weighted by defining a viscosity-related length d^ = 



(//f/of 1 AP~ 1 ) 1 / 2 and flow rate Q^ = (/ifp^AP -1 ) 1 / 2 . Using these and equation (10), an 
entirely analogous expression is obtained in the limit of dominant inertia, i. e. when dj » 
dy, and Qi » Q M : 

dj/dy = UQt/Q^ 2 . (12) 

This equation expresses the jet diameter as compared to a limit where viscous effects become 
important. Again, jet diameters and flow rates comparable to d^ and Q^, respectively, 



amount to non-negligible viscous effects and significant departures from predictions (11) or 



(J12J). 

A third expression can be obtained in terms of the orifice diameter D, and the maximum 
liquid flow rate that can be ejected through the orifice for a given AP in the absence of 
viscous effects: Q ma x = ( 7 >" 2 /8) 1 ' /2 Q m , where -naturally- Q m = (D^APpJ 1 ) 1 / 2 is obtained 
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from dimensional analysis using {AP,D,a}. Using (10) anew, one has: 



dj/D = (Qi/Q max ) 1/2 = k d (Qi/Q m ) 1/2 . (13) 

This alternative expression reflects how close the experiment is from a situation where the 
entire orifice section is filled with liquid: it provides information -from continuity arguments- 
on the fraction of the orifice cross section occupied by the liquid jet. 



Each of the above three expressions (11), (12) and (13) amount to interesting but partial 



pictures of the particular working conditions of our system in a given flow situation. Taken 
as a whole, they provide a more complete picture on the FF jetting conditions. Some 
corrections can be obtained for several neglected effects. [2T] 

1. Correction for surface tension effects 

The liquid surface tension reduces the effective pressure drop APi in the liquid stream as 

APi = AP - 2a/ dj. (14) 

Consequently, the jet velocity decreases and its diameter increases accordingly. The re- 
sulting expression for the non-dimensional jet diameter dj/d , neglecting third order terms 
proportional to 0(d /dj) « 1, reads: 

dj/do = (8/7r 2 ) 1/4 (^/Qo) 1/2 + 1/2. (15) 

In other words, the second order correction of the jet diameter dj to account for surface 
tension effects is asymptotically equal to d /2. 

2. Correction for liquid viscosity effects (extensional stresses) 

Assuming that the extensional viscous forces in the liquid are smaller than inertia, the 
balance of the different terms of the momentum equation, including the second order terms 
of the expansion, leads to the following order of magnitude for the correction to the first 



order diameter (10): 



d P = 



1/2' 

d. 



^iraax 



M V Qi 



(16) 
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3. Correction for tangential stresses owing to the gas stream 

In the same way, the diameter correction (decrease) owing to the momentum injected by 
the much faster gas stream through the jet surface is of the order of: 

«. = °(^f (it) 

where p g and U g are the gas viscosity and velocity. The latter is of the order of U g ~ 

1 /9 

O (AP/p g ) ' , where p g is the gas density 

The relative weight of these three corrections provides information on the importance of 
the surface tension and the viscosity of the liquid and gas phases. Interestingly, for most 
common solvents, these relative weights are of order unity. This happens to be the case 
when measuring the relative importance of the surface tension and the gas tangential stress 
effects for water focused by any gas at standard conditions. Therefore, since both corrections 
are opposite, the best agreement with experimentally measured jet diameters and numerical 



simulations is obtained, interestingly enough, using the first order expression (10), or its 



alternative forms (11 13) 



4- Correction owing to the nozzle flow pattern 

The jet diameter as measured at the nozzle may also differ from the simplest theoretical 



prediction given by Eq. ( 10 ) because of local flow effects. A complex but symmetric structure 
develops owing to the coexistence of (i) a core potential flow and (ii) the detachment of a 
radially convergent boundary layer at the inner lip of the nozzle. In any real situation 
where the gas viscosity is non-zero and the continuum hypothesis holds, this flow pattern 
is not aptly described by the pure potential flow through a round orifice given by Morse 
and Feshback[55j (page 1294) for a stationary discharge. The potential flow solution is 
characterized by an axial velocity distribution with a minimum value at the axis, v(r = 0) = 
2Q g /(7cD 2 ) (half the average velocity through the orifice), and an infinite value at r = D/2, 
Q g being the theoretical gas flow rate discharged. The actual flow geometry is characterized 
by the well known vena contracta effect, a consequence of the radial momentum carried 
by the collapsing potential flow, which slips at the nozzle border owing to the boundary 
layer. The vena contracta flow exhibits an axial velocity distribution which echoes the 



31 



potential flow solution, showing a local minimum velocity at the axis, and a maximum 
value at the streamlines coming just from the outside of the boundary layer detached at the 
orifice (see Fig. [9]) . The immediate consequence of this particular flow structure is that the 
transversal pressure gradients are negligible only sufficiently far downstream of the inner 
lip of the nozzle: in fact, they become negligible at the axial downstream station where 
the vena contracta effect ends, i.e. where the streamlines become almost parallel. It occurs 
relatively close to the inner orifice plane, at an approximate D/2 downstream distance. 
From this point downstream (before shear instabilities of the gas stream with the external 
environment develop), the gas pressure can be considered almost constant, equal to the 
outside stagnation value. It is at this point where the liquid jet diameter obtained from the 



numerical simulation should be compared to the simplest prediction (10). 



A. Scaling of the recirculation zone 

For a given gas flow rate Q g and orifice diameter D = 2R, the typical gas velocity close to 
the meniscus surface can be estimated as V = Q g /{j{R 2 ). Given the small p = p g /pi values 
in liquid jets focused by gas, liquid velocities are much smaller than V everywhere. As the 



liquid approaches the neck, the boundary layer will collapse (Fig. 22). This implies that at 
least a liquid flow rate 

Qr ~ U s 5 2 (18) 

would be drawn into the jet in the absence of recirculation (U s is the velocity of the interface, 
that can be obtained from V, and the densities and the viscosities ratios[55]). On the 
contrary, whenever Q L < Q R , part of Qr must have been recirculated back into the meniscus 



(Fig. 22). Therefore Qr can be interpreted as the minimum flow rate for no recirculation 
(scaled as Q r = Qr/Q 9 ). 

The boundary layer in the liquid meniscus is confined. It grows along the cone during 
lengths comparable to R (the orifice radius) till the apex of the meniscus is reached. In this 
area, the gas speed gradients are steep: 

5 t ~ {pi l R/p l U s ) 1 l 2 . (19) 

Whenever there is recirculation, the peripheral boundary layers merge at the meniscus apex 
and give rise to a jet, whose initial radius at the neck will accordingly be of the same order. 
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D = 2R 



FIG. 22: Sketch of the recirculation zone, showing boundary layers, cell size (S r ), and typical 
velocities. 

In the absence of liquid emission, maximum recirculation will be observed. Experimentally, 
however, a dripping instability will occur before reaching this limit. In the opposite case (no 
recirculation), the boundary layers do not merge, and an inviscid core should be observed at 
the neck. The threshold flow rate for recirculation can therefore be estimated as Qr ~ U s 5f, 
a result which happens to be independent of the gas velocity. In effect, by definition of 
the meniscus boundary layer, the viscous stress niU s /5f must be of the same order as the 
momentum convection piU^/R, so that, interestingly: 

,M . „ P 



Qi 



R- 



Wr 



(20) 



Pi I^Re 

This scaling is fully confirmed by the numerical simulations: the values of Q r ■ Re are 0.6768 
for case 1 and 0.6596 for case 2, deviating by less than 2.6% from the scaling predictions. 

Assume now the recirculation cell to be S r in axial length. The backflow Qb = Qr — Qi 
will come to rest within a length of the order S r . In this length, viscous momentum diffusion 
should slow down the flow and deflect both the incoming flow injected by the feeding tube 
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and the recirculated flow at the axis (Fig. 22). Thus, viscous and inertia forces should 



balance within that length S r : in other words, the liquid Reynolds number associated to 
axial lengths of order S r should be of order unity so as to deflect the unidirectional flow 
issuing from the feed tube (Hagen-Poiseuille). This is in analogy to the entry length or exit 
length in laminar pipe flow. Two cases need to be considered, depending of the relative size 
of the cell compared to the feed tube radius Ri. 

• When S r < R\, viscous stress, of the order 0(piQbS^ 3 ), balances inertia, O^piQ^S' 4 "), 
which leads to S r ~ PiQbPi^ '■ 

• When S r > R\, viscous stress, 0(piQbRi 4 ), balances inertia, O^piQgR^S' 1 ), leading 
again to S r ~ PiQbUi "■ 

Interestingly enough, again, the length of the recirculation flow is independent of the gas 
flow for any given geometry. The latter scaling can be expressed in non-dimensional terms 
as: 

Sr^Sr/R-ptiQR-Q^p^R- 1 . (21) 



Using equation (18), one may write: 

s r = C l - C 2 Re Rl (22) 

where Re R = piQif(piR) is a Reynolds number of the liquid flow, and C\ and C 2 are 

constants which depend on the geometry only (i.e., R\/R, H/R, etc.). In our case, we have 



represented all our measured s r values from numerical simulations versus Rcr in Fig. ( 23 ) 



Linear fitting to all points leads to C\ = 2.636 and C 2 = 0.0819 with a correlation coefficient 



of 94.4%. Equation (22) can be expressed in terms of Q r — Q and Re as well, as: 

s r = kpp-\Q r -Q), (23) 

where k is again a constant depending on the geometry only, in full agreement with expression 
(|8J), as anticipated by experiments. 



VI. EXPERIMENTS 

In the following, we provide experiments corresponding to the same local geometrical 
parameters as in cases 1 and 2 in the vicinity of the exit orifice. The basic flow focusing 
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FIG. 23: Recirculation cell size s r as a function of Rcr: dots, squares and stars are obtained by 
numerical simulation; the line is a theoretical prediction resulting from dimensional arguments. 
An additional series of simulations have been performed for an intermediate gas flow condition 
(Re = 698.75 and We = 18.31) to assess the validity of the scaling proposed: note the good 
degree of collapse obtained. The small deviations can be attributable to the small differences in 
the geometry of the cone for different gas flow conditions. 

chamber is a box consisting of five aluminum faces and one clear methacrylate face. It is 5 
cm by 5 cm by 5.65 cm, with its longest side along the capillary/orifice axis. The chamber 
is situated with the methacrylate face horizontal and pointing upwards, the capillary being 
located parallel to this face. The orifice is made in a stainless steel orifice disk attached to 
the box side, perpendicular and opposite to the capillary. The disk is 4.0 mm in diameter 
with a thickness of 75 /xm and an orifice of diameter 0.200 mm. Both the air tube and 
the capillary enter through the face opposite the orifice. After the capillary tube has been 
aligned with the orifice, the distance H from the tube to the orifice can be simply adjusted 
by carefully sliding the capillary in its housing on the opposite face to the orifice disk. H 



is measured with a microscope through the methacrylate face. Fig. 24 shows some views 
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FIG. 24: (Left) Experimental tube-orifice set up as numerically simulated in this work {D = 
200/nn, Di = 150/xm, H = lOO^m; here, AP = lOiTPa, Qi = 3 mL/h). (Right) Photographs 
of experimental conditions with twice the distance from the feeding tube to the exit orifice, using 
a different tube material (fused silica): (A) jetting {D = 200/im, D\ = 150/um, H = 200/nn, 
AP = 30KPa, Qi = Q.l mL/h) and (B) dripping (as in A, with Qi = 2.8 mL/h). 



of the feeding tube-orifice setup as seen through the thick methacrylate window (inevitable 
liquid spills leave behind some debris on the inner face of the window causing a blurred 



image). In particular, Fig. 24 (left) shows the geometry numerically simulated in this work. 



After setting H and ensuring that the capillary is perfectly coaxial with the nozzle orifice, 
the pressure is set using a pressure gauge and a pressure meter. A water flow rate is then 
supplied using a syringe pump (Cole-Palmer 74900 Series) with a 20 ml syringe. The system 
is given sufficient time to relax until either a characteristic steady or unsteady flow is present. 
This can be checked by illuminating the jet that exits the orifice or by looking at the meniscus 
when the distance H is 0.100 mm or greater. Unsteady jet flow appears very faint to the 
naked eye and contains thin streaks of water along with large scattered spray. This is in 
significant contrast to steady jet flow, which has bright illumination as a result of a finer, 
concentrated stream with uniform characteristics. In experiments where the meniscus was 



visible, it was also possible to discriminate steady versus unsteady flow (see Fig, 24 (right) 



A -jetting- and B -dripping-), in perfect correlation with the spray observations: a steady 



meniscus had sharp edges and a clear, unwavering glass-like appearance (Fig. 24 A, see 



steady jet reflected in the metal plate), while an unsteady meniscus had blurred edges and 



flickered (Fig. 24 B, no jet is visible at all). Both the jet test and meniscus test displayed 
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FIG. 25: Jetting -dripping transition in the {Rei,Wei} plane. Diamonds: experimentally deter- 
mined conditions (filled symbols: liquid flow rate decreasing -"Down"; open symbols: liquid flow 
rate increasing- "Up" . In most cases, both "Up" and "Down" points coincide) . Circles: numerically 
tested conditions. Filled circles: jetting conditions. Open circles: dripping conditions. 

clear and abrupt transitions between the two states. Once unsteady flow is established 
for a given pressure, the rate determined by the syringe pump is increased in steps of 0.1 
ml/hr. After each flow rate increase, a 30 seconds waiting period was established, so as to 
ensure that the system had relaxed and all the readings were accurate. This period has been 
chosen after it was found that 15s was not enough to observe fluctuations in the system: 
occasionally a steady regime would revert back to an unsteady one after the 15s period. The 
30s delay has proven long enough to accurately characterize the flow. Accordingly, the rate 
was increased until the unsteady jet sharply transitioned to a steady one; at this point, the 
flow rate was read from the syringe pump and recorded as the minimal flow rate (increasing, 



or "up"; Fig. 25). Keeping this same steady flow rate the process is then reversed to find 



the minimum decreasing (or "down") flow rate (steps of 0.1 ml/hr and intervals of 30 s 
until an unsteady regime developed). When the flow became unsteady, a rate 0.1 ml/hr 
above the reading on the syringe pump was recorded, since the rate which produced the last 



37 



steady flow (i.e. minimum flow rate) was one step (0.1 ml/hr) higher. The resulting value 
was recorded as the minimal (decreasing, or "down" ) flow rate. This process is repeated for 
varying pressures and distances of H to get an accurate mapping of minimal jetting flow 
rates as a function of varying geometry and flow conditions. Following this procedure, we 



collected the experimental data plotted in Fig. 25 for H/R = 1. The gas (air) pressure AP 
increases as indicated by the arrow. 

Six conditions numerically tested for cases 1 and 2 are plotted in Figj25j In order to 
make our results readily translatable in most of the capillary jet stability literature (which 
uses the jet radius as a characteristic length), we may introduce liquid Reynolds and Weber 



numbers consistent with previous definitions and using scaling law (10): 



As it follows from the plot, using these definitions, jetting or dripping conditions are accu- 
rately predicted by the numerical model. This lends additional support to the use of full 
VOF simulation analysis in flow focusing systems. 

VII. CONCLUSIONS 

The cone-jet geometry associated with flow focusing has been handled by a diversity 
of tools, numerical, experimental and theoretical. Order-of-magnitude estimations follow 
from dimensional arguments: such procedures contribute a valuable theoretical framing and 
provide the scaling criteria for data representation. Analytical approaches are generally 
based on the consideration of a perfectly cylindrical infinite jet, a simplification that ignores 
the influence of the meniscus (a source of instability) and the role of streamline convergence 
or divergence in the jet. Experiments are burdened by the diversity of influencing parameters 
and visualization difficulties associated with the small scale of the meniscus and jet. 

In this paper, experimental results are backed up by a numerical simulation based on VOF 
elements. Numerical schemes allow a more systematic exploration of the parametric influ- 
ence. In addition, the shortcomings of theoretical models (unavoidable in a situation where 
the geometry of the fluid domain is complex, as in a cone-jet flow pattern) are overcome, 
and a detailed description of the streamlines can be readily obtained. 

The key results of the above exploration are the following: 
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The theoretical scaling leading to jet diameter estimates is confirmed by the simulation. 
The expressions for flow focusing scales, notwithstanding their simplicity, are therefore 
to be considered a reliable shortcut for the prediction of jet dimensions. 

The complete sequence from meniscus growth to jet emission (jetting regime) and to 
the sequential filling of drops (dripping regime) is portrayed in detail. 

The jetting-dripping transition is documented in detail, both by experiment and sim- 
ulation. A two-branch structure is observed in the plot, showing the simultaneous 
influence of the jet and the meniscus as instability sources. Incipient dripping (Fig. 
fi~3|) is shown to give rise to highly irregular fluctuations; while fully developed dripping 



(Fig. 15) produces perfect cycles of drop detachment. 



• A recirculation cell is identified in the jetting regime at the meniscus tip. This oc- 
currence appears to be linked to intensive forcing by the gas sheath, leading to high 
interface velocity along the meridians; the issuing jet is unable to convey all of the 
mobilized flow, so that a return flow around the axis is observed. The recircula- 
tion cell grows as the liquid flow rate is reduced: eventually, dripping conditions are 
reached. Similar recirculation cells have been observed in electrospray cones, under 
thread emission, and in liquid-liquid two-dimensional flow focusing, assisted or not by 
a surfactant [SJ [2HJ EH]. All of the recirculation instances reported thus far appear 
to share a common attribute: strong interfacial forcing, either electric, capillary or 
hydrodynamic. 

• A reliable scaling is provided, identifying the parametric conditions where recirculation 
is to be expected and estimating the size and flow rates of the cell. 

A key feature in the flow pattern explored is the recirculation cell, and its conceptual 
link to the merging of the boundary layers which grow from the meniscus edge and fuse 
together at the neck of the jet. Controllable recirculation is an extremely attractive feature, 
providing adjustable residence times within a very simple flow setup. The cell can be viewed 
as a flow trap or reactor, where biosynthesis or chemical operations take place in a protected 
environment; the liquid flow rate can be increased to flush the recirculation products. 

An additional focus deals with the peculiarities of the jetting and dripping regimes under 
the influence of a co-flowing sheath current. The aspect of the jet and droplet train and the 
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dynamics of the meniscus (an indicator of dripping) are a contribution to a problem whose 
complexity forbids a global theoretical approach. 
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